double  interaccion(caja *bx,double rc)
{
	int ind;
	int ank, ank2, ankh;
	double virial=0.0;
	ank=6;
	ank2=2*ank;
	ankh=ank/2;
	double dist2,rqInv,dr[3], tsn, ccel3,potdm, ePoten;
	double conrep= 4.0*kB*pow(1.0,ank2); //Parte repulsiva del potencial Leonard-Jones
	double conref= ank2*conrep; //Parte repulsiva de la fuerza del Leonard-Jones
	double conat= 4.0*kB*pow(1.0,ank); //Parte atractiva del potencial Leonard-Jones
	double conatf= ank*conat;// Parte atractiva de la fuerza del Leonard-Jones

	//Ponemos las fuerzas de las particulas de la caja a 0
	bx->initForces();

	//Correccion al cutoff del potencial
	double LJcut= pow(1.0/(rc*rc),ankh)*(pow(1.0/(rc*rc),ankh)*conrep-conat);

	ePoten=0.;

	for(int i = 0;i<bx->npart-1;i++)
	{
		//Recorremos los vecinos de la particula i
		for(int j=0;j<bx->part[i].nvecinos;j++)
		{
			ind = bx->part[i].lstVecinos[j];
			dist2= bx->distancia2(i,ind);
			rqInv =1.0/dist2;//1/dist^2
			tsn=pow(rqInv,ankh);
			//Calculo fuerza ccel3=F/r=-(dPOTDM/dr)/r 
			ccel3 =rqInv*tsn*(10.0*conref*tsn-conatf);
			potdm = tsn*(conrep*tsn-conat);
			//Actualizamos la energia potencial
			ePoten+=potdm;
			//Calculamos los componentes de las fuerzas
			for (int k=0;k<3;k++)
			{
				dr[k]= bx->dr[k];
				//Para la particula i
				bx->part[i].force[k]+=ccel3*dr[k];
				//Para la particula ind
				bx->part[ind].force[k]-=ccel3*dr[k];
			}
			//Actualizamos el virial
			 virial+=ccel3*(-dr[0]*dr[0]-dr[1]*dr[1]-dr[2]*dr[2]);
		}
	}
	//Introducimos la correccion del cutoff del potencial
	ePoten-=LJcut*(bx->npart)*((bx->npart)-1)/2.0;
	bx->epoten=ePoten;
	bx->virial=virial;
	return ePoten;

}

double interaccionLinkCellList(caja *bx)
{
	int c, c1;
	int iii,jjj;
	int shift[3];
	int ank, ank2, ankh;
	int count=0;//Numero de pares interactuantes
	ank=6;
	ank2=2*ank;
	ankh=ank/2;
	double dist2,rqInv,dr[3], tsn, ccel3,potdm, ePoten;
	//double conrep= 4.0*pow(1.0,ank2); //Parte repulsiva del potencial Leonard-Jones
	//double conref= ank2*conrep; //Parte repulsiva de la fuerza del Leonard-Jones
	//double conat= 4.0*pow(1.0,ank); //Parte atractiva del potencial Leonard-Jones
	//double conatf= ank*conat;// Parte atractiva de la fuerza del Leonard-Jones
	double conrep= 4.0; //Parte repulsiva del potencial Leonard-Jones
	double conref= 48.0; //Parte repulsiva de la fuerza del Leonard-Jones
	double conat= 4.0; //Parte atractiva del potencial Leonard-Jones
	double conatf= 24.0;// Parte atractiva de la fuerza del Leonard-Jones

	double rc = bx->rcut;
	double rc2=rc*rc;
	double virial=0;
	//Ponemos las fuerzas de las particulas de la caja a 0
	bx->initForces();
	//Correccion al cutoff del potencial
	double LJcut= pow(1.0/(rc*rc),ankh)*(pow(1.0/(rc*rc),ankh)*conrep-conat);

	ePoten=0.;
	//Recorremos todas las celdas
	for (int ix=0;ix<bx->ncell;ix++)
	{
		for (int iy=0;iy<bx->ncell;iy++)
		{
			for (int iz=0;iz<bx->ncell;iz++)
			{
				///Calculamos el indice de la celda (ix,iy,iz)
				c = int(ix*bx->ncell2+iy*bx->ncell+iz);
				///Recorremos las particulas de esta celda y las celdas vecinas
				for (int ixred=ix-1;ixred<=ix+1;ixred++)
				{
					/*shift: condiciones de contorno periodicas tenidas en cuenta para
					localizar las celdas vecinas*/
					if (ixred < 0)
					{
						shift[0] = int(ixred+bx->ncell);
					}
					else
					{
						if (ixred>=bx->ncell)
						{
							shift[0] = int(ixred-bx->ncell); 
						}
						else
						{
							shift[0] = int(ixred);
						}
					}
					for (int iyred=iy-1;iyred<=iy+1;iyred++)
					{
						if (iyred < 0)
						{
							shift[1] = int(iyred+bx->ncell);
						}
						else
						{
							if(iyred>=bx->ncell)
							{
								shift[1] = int(iyred-bx->ncell);
							}
							else
							{
								shift[1] = int(iyred);
							}
						}
						for (int izred=iz-1;izred<=iz+1;izred++)
						{
							if (izred < 0)
							{
								shift[2] = int(izred+bx->ncell);
							}
							else
							{
								if(izred>=bx->ncell)
								{
									shift[2] = int(izred-bx->ncell);
								}
								else
								{
									shift[2] = int(izred) ;
								}
							}
							///Calculamos el indice c1 de la celda vecina
							c1 = int(((shift[0]+bx->ncell)%bx->ncell)*bx->ncell2+((shift[1]+bx->ncell)%bx->ncell)*bx->ncell
								+((shift[2]+bx->ncell)%bx->ncell));
							///Buscamos la primera particula iii de la celda c
							iii = bx->head[c];
							while (iii != -1) ///iii=-1 ->No hay particulas en la celda
							{
								///Buscamos la primera particula jjj de la celda c1
								jjj =bx->head[c1];
								while (jjj != -1) ///jjj=-1->No hay particulas en la celda
								{
									//Evitamos el conteo doble de las interacciones
									if (iii < jjj) 
									{
										///Calculamos la distancia con condiciones periodicas entre iii y jjj
										dist2= bx->distancia2(iii,jjj);
										//Si la distancia es menor que el cutoff calculamos la interaccion
										if (dist2 < rc2)
										{
											count++;
											///Calculamos la fuerza entre iii y jjj
											rqInv =1.0/dist2;//1/dist^2
											tsn=pow(rqInv,ankh);//1/dist^6
											//ccel3 = F/r=-(dPOTDM/dr)/r 
											ccel3 =rqInv*tsn*(conref*tsn-conatf);
											potdm = tsn*(conrep*tsn-conat);
											//Acualizamos la energia potencial
											ePoten+=potdm;
											//Calculamos los componentes de las fuerzas
											for (int k=0;k<3;k++)
											{
												dr[k]= bx->dr[k];
												//Para la particula iii
												bx->part[iii].force[k]+=ccel3*dr[k];
												//Para la particula jjj
												bx->part[jjj].force[k]-=ccel3*dr[k];
											}
											//Actualizamos el virial
											virial-=ccel3*(dr[0]*dr[0]-dr[1]*dr[1]-dr[2]*dr[2]);
										}
									}
									//Pasamos a la siguiente particula de la celda c1
									jjj = bx->link[jjj];
								}
								//Pasamos a la siguiente particula de la celda c
								iii = bx->link[iii];
							}
						}
					}
				}
			}
		}
	}
	//Introducimos la correccion del cutoff del potencial
	ePoten-=LJcut*count;
	bx->epoten=ePoten/bx->npart;
	bx->virial=virial;
	return ePoten;
}


double interaccionLinkCellList(int ind, caja *bx)
{
	int c, c1;
	int iii,jjj;
	int shift[3];
	int ank, ank2, ankh;
	ank=6;
	ank2=2*ank;
	ankh=ank/2;
	double dist2,rqInv,dr[3], tsn, ccel3,potdm, ePoten;
	double conrep= 4.0*pow(1.0,ank2); //Parte repulsiva del potencial Leonard-Jones
	double conref= ank2*conrep; //Parte repulsiva de la fuerza del Leonard-Jones
	double conat= 4.0*pow(1.0,ank); //Parte atractiva del potencial Leonard-Jones
	double conatf= ank*conat;// Parte atractiva de la fuerza del Leonard-Jones
	double rc = bx->rcut;
	double rc2=rc*rc;
	double virial=0;
	//Ponemos las fuerzas de las particulas de la caja a 0
	bx->initForces();
	//Correccion al cutoff del potencial
	double LJcut= pow(1.0/(rc*rc),ankh)*(pow(1.0/(rc*rc),ankh)*conrep-conat);

	ePoten=0.;
	//Recorremos todas las celdas
	for (int ix=0;ix<bx->ncell;ix++)
	{
		for (int iy=0;iy<bx->ncell;iy++)
		{
			for (int iz=0;iz<bx->ncell;iz++)
			{
				///Calculamos el indice de la celda (ix,iy,iz)
				c = int(ix*bx->ncell2+iy*bx->ncell+iz);
				///Recorremos las particulas de esta celda y las celdas vecinas
				for (int ixred=ix-1;ixred<=ix+1;ixred++)
				{
					/*shift: condiciones de contorno periodicas tenidas en cuenta para
					localizar las celdas vecinas*/
					if (ixred < 0)
					{
						shift[0] = int(ixred+bx->ncell);
					}
					else
					{
						if (ixred>=bx->ncell)
						{
							shift[0] = int(ixred-bx->ncell); 
						}
						else
						{
							shift[0] = int(ixred);
						}
					}
					for (int iyred=iy-1;iyred<=iy+1;iyred++)
					{
						if (iyred < 0)
						{
							shift[1] = int(iyred+bx->ncell);
						}
						else
						{
							if(iyred>=bx->ncell)
							{
								shift[1] = int(iyred-bx->ncell);
							}
							else
							{
								shift[1] = int(iyred);
							}
						}
						for (int izred=iz-1;izred<=iz+1;izred++)
						{
							if (izred < 0)
							{
								shift[2] = int(izred+bx->ncell);
							}
							else
							{
								if(izred>=bx->ncell)
								{
									shift[2] = int(izred-bx->ncell);
								}
								else
								{
									shift[2] = int(izred) ;
								}
							}
							///Calculamos el indice c1 de la celda vecina
							c1 = int(((shift[0]+bx->ncell)%bx->ncell)*bx->ncell2+((shift[1]+bx->ncell)%bx->ncell)*bx->ncell
								+((shift[2]+bx->ncell)%bx->ncell));
							///Buscamos la primera particula iii de la celda c
							iii = bx->head[c];
							while (iii != -1) ///iii=-1 ->No hay particulas en la celda
							{
								///Buscamos la primera particula jjj de la celda c1
								jjj =bx->head[c1];
								while (jjj != -1) ///jjj=-1->No hay particulas en la celda
								{
									//Evitamos el conteo doble de las interacciones
									if (iii < jjj) 
									{
										///Calculamos la distancia con condiciones periodicas entre iii y jjj
										dist2= bx->distancia2(iii,jjj);
										//Si la distancia es menor que el cutoff calculamos la interaccion
										if (dist2 < rc2)
										{
											///Calculamos la fuerza entre iii y jjj
											rqInv =1.0/dist2;//1/dist^2
											tsn=pow(rqInv,ankh);
											//ccel3 = F/r=-(dPOTDM/dr)/r 
											ccel3 =rqInv*tsn*(conref*tsn-conatf);
											potdm = tsn*(conrep*tsn-conat);
											//Acualizamos la energia potencial
											ePoten+=potdm;
											//Calculamos los componentes de las fuerzas
											for (int k=0;k<3;k++)
											{
												dr[k]= bx->dr[k];
												//Para la particula iii
												bx->part[iii].force[k]+=ccel3*dr[k];
												//Para la particula jjj
												bx->part[jjj].force[k]-=ccel3*dr[k];
											}
											//Actualizamos el virial
											virial+=ccel3*(dr[0]*dr[0]-dr[1]*dr[1]-dr[2]*dr[2]);
										}
									}
									//Pasamos a la siguiente particula de la celda c1
									jjj = bx->link[jjj];
								}
								//Pasamos a la siguiente particula de la celda c
								iii = bx->link[iii];
							}
						}
					}
				}
			}
		}
	}
	//Introducimos la correccion del cutoff del potencial
	ePoten-=LJcut*(bx->npart)*((bx->npart)-1)/2.0;
	bx->epoten=ePoten;
	bx->virial=virial;
	return ePoten;
}
